Celem raportu było określenie głównych przyczyn stopniowego spadku rozmiaru śledzia oceanicznego wyławianego w Europie na przestrzeni ostatnich lat. Dostarczony zbiór danych zapewniał informacje o następujących parametrach:
W ramach połowu jednej jednostki wybierano od 50 do 100 sztuk trzyletnich śledzi. Analiza wykazała, że kluczowe znaczenie dla długości śledzi mają temperatura przy powierzchni wody, natężenie połowów w regionie oraz wzrost oscylacje północnoatlantyckie. Wzrost oscylacji powodował malenie długości śledzia. Zmniejszenie się długości śledzi mogło być też spowodowane niedoborem planktonu gatunków chel1 i lcop1. Zwiększanie natężenia połowów w regionie powodowało zwiększanie długości śledzi - prawdopodobnie dlatego, że zwiększała się ilość pokarmu przypadającego na pozostałe w wodzie jednostki.
Analizy dokonano w środowisku RStudio z wykorzystaniem dokumentu R Markdown i z użyciem paczki knitr. Skorzystano też z następujących bibliotek:
Powtarzalność wyników dla operacji losowych zapewnia ustawienie stałego ziarna.
set.seed(1)
Dane wczytane z pliku sledzie.csv będą w dalszej kolejnosci przetwarzane w postaci typu data frame. Podczas wczytywania zbioru konieczne jest określenie typu danych dla każdej kolumny jako numeryczny.
fname = "sledzie.csv"
headset <- read.csv(fname, header = TRUE, nrows = 10)
dataset <- read.csv(fname, header = TRUE, na.strings = "?", colClasses = "numeric")
df <- as.data.frame(dataset)
sapply(df, class)
## X length cfin1 cfin2 chel1 chel2 lcop1
## "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## lcop2 fbar recr cumf totaln sst sal
## "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## xmonth nao
## "numeric" "numeric"
Fragment zbioru danych przedstawia się następująco:
head(df)
## X length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr
## 1 0 23.0 0.02778 0.27785 2.46875 NA 2.54787 26.35881 0.356 482831
## 2 1 22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 3 2 25.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 4 3 25.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 5 4 24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 6 5 22.0 0.02778 0.27785 2.46875 21.43548 2.54787 NA 0.356 482831
## cumf totaln sst sal xmonth nao
## 1 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 2 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 3 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 4 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 5 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 6 0.3059879 267380.8 14.30693 35.51234 7 2.8
Dane w kolumnach nie zmieniają się z dużą dynamiką, dlatego wiersze z brakującymi informacjami nie zostaną pominięte, a uzupełnione, gdyż ta operacja na pewno nie wprowadzi do zbioru dużych perturbacji. Do uzupełnienia pominiętych wartości użyta zostanie wartość średnia kolumn. Nie jest to najlepsze możliwe rozwiązanie, ale jego znaczną zaletą jest szybkość działania. W tym przypadku średnia jest metodą akceptowalną z uwagi na przypuszczalnie niską wariancję zmiennych.
sapply(df, function(x) sum(is.na(x)))
## X length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr
## 0 0 1581 1536 1555 1556 1653 1591 0 0
## cumf totaln sst sal xmonth nao
## 0 0 1584 0 0 0
sum(complete.cases(df))
## [1] 42488
Zbiór zawiera 42488 kompletnych przypadków pośród 52582 rekordów. Brakujące wartości występują w kolumnach cfin1, cfin2, chel1, chel2, lcop1, lcop2 oraz sst, a ich liczba jest do siebie zbliżona. Poniższy wykres prezentuje rozkład brakujących zmiennych.
missing_plot(df)
Braki informacji są przypadkowe, często dla danego przypadku brakuje tylko jednej zmiennej.
clean_df <- na_mean(df)
anyNA(clean_df)
## [1] FALSE
Operacja prawidłowo uzupełnia zbiór, nie pozostawiając brakujących wartości.
Zbiór zawiera 52582 rekordy.
Podstawowe statystyki na temat wartości kolumn:
summary(clean_df)
## X length cfin1 cfin2
## Min. : 0 Min. :19.0 Min. : 0.0000 Min. : 0.0000
## 1st Qu.:13145 1st Qu.:24.0 1st Qu.: 0.0000 1st Qu.: 0.2778
## Median :26291 Median :25.5 Median : 0.1333 Median : 0.7012
## Mean :26291 Mean :25.3 Mean : 0.4458 Mean : 2.0248
## 3rd Qu.:39436 3rd Qu.:26.5 3rd Qu.: 0.3603 3rd Qu.: 1.9973
## Max. :52581 Max. :32.5 Max. :37.6667 Max. :19.3958
## chel1 chel2 lcop1 lcop2
## Min. : 0.000 Min. : 5.238 Min. : 0.3074 Min. : 7.849
## 1st Qu.: 2.469 1st Qu.:13.589 1st Qu.: 2.5479 1st Qu.:17.808
## Median : 6.083 Median :21.435 Median : 7.1229 Median :25.338
## Mean :10.006 Mean :21.221 Mean : 12.8108 Mean :28.419
## 3rd Qu.:11.500 3rd Qu.:27.193 3rd Qu.: 21.2315 3rd Qu.:37.232
## Max. :75.000 Max. :57.706 Max. :115.5833 Max. :68.736
## fbar recr cumf totaln
## Min. :0.0680 Min. : 140515 Min. :0.06833 Min. : 144137
## 1st Qu.:0.2270 1st Qu.: 360061 1st Qu.:0.14809 1st Qu.: 306068
## Median :0.3320 Median : 421391 Median :0.23191 Median : 539558
## Mean :0.3304 Mean : 520367 Mean :0.22981 Mean : 514973
## 3rd Qu.:0.4560 3rd Qu.: 724151 3rd Qu.:0.29803 3rd Qu.: 730351
## Max. :0.8490 Max. :1565890 Max. :0.39801 Max. :1015595
## sst sal xmonth nao
## Min. :12.77 Min. :35.40 Min. : 1.000 Min. :-4.89000
## 1st Qu.:13.63 1st Qu.:35.51 1st Qu.: 5.000 1st Qu.:-1.89000
## Median :13.86 Median :35.51 Median : 8.000 Median : 0.20000
## Mean :13.87 Mean :35.51 Mean : 7.258 Mean :-0.09236
## 3rd Qu.:14.16 3rd Qu.:35.52 3rd Qu.: 9.000 3rd Qu.: 1.63000
## Max. :14.73 Max. :35.61 Max. :12.000 Max. : 5.08000
Widzimy, że długość śledzi waha się w przedziale [19.0, 32.5] i ma średnią 25.3 oraz zbliżoną do niej medianę - 25.5. Statystyki dostępności planktonu wykazują mniejszą medianę od średniej (z wyjątkiem planktonu chel2). Znaczy to, że w pewnym okresie plankton dostępny był w bardzo dużej ilości, ale przez większość czasu jego ilość była znacznie bardziej ograniczona. Największą różnicę widać w przypadku planktonu cfin2. Temperatura przy powierzchni wody (sst) oraz zasolenie oceanu (sal) prezentują małe wahania wartości, jednak dalsze analiza dowiedzie, że to właśnie te parametry mają krytyczne znacznie dla długości śledzi.
W pierwszej kolejności przeanalizowane zostaną rozkłady zmiennych z użyciem wykresów typu Box plot i histogramów. Grafy rozkładów zostały pogrupowane pod kątem zakresu wartości. Kod generujący wykresy, z powodu długości, nie zostanie zaprezentowany w raporcie.
Na powyższym wykresie możemy zauważyć przede wszystkim, że większe długości śledzi występowały dla większych ilości planktonu chel1 i lcop1. Dla większych długości śledzi malejącą tendencję miały też oscylacje północnoatlantyckie. Wygląda też na to, że miesiąc połowu nie miał prawie żadnego wpływu na długość śledzia. Najmniejszej długości śledzi występowały, gdy dostępność planktonu wszystkich gatunków była minimalna.
Możemy zauważyć, że dla większości zmiennych mediana nie byłaby dobrą metodą zastępowania brakujących wartości, gdyż zwykle jest przesunięta ku któremuś z krańców rozkładu. W zbiorze często pojawiają się wartości zmiennych mniejsze lub większe od mediany. W ogólności występują pojedyncze obserwacje odstające, co znaczy, że zbiór nie jest w znacznym stopniu zaszumiony. W oceanie dominowała dostępnosć planktonu lcop2, który jednocześnie miał największe wahania wartości, ale jego mediana nie była w znacznym stopniu przesunięta ku jednemu z krańców rozkładu. Drugim dominującym źródłem pokarmu był plankton chel2. Plankton lcop1 miał duże rozkłady wartości, ale zwykle było go niewiele - mediana jego jest zbliżona do mediany planktonu chel1.
Na histogramach widać jeszcze wyraźniej, że przez większość czasu dostepność planktonu wszystkich rodzajów była stosunkowo mała w porównaniu do najwyższych zanotowanych wartości. Największy niedobór planktonu został zanotowany dla gatunków cfin1, lcop2 i lcop1. Połowów zwykle dokonywano w miesiącach letnich (od lipca do października). Zasolenie w przeważającej mierze utrzymywało się na średniej wartości, jednakże zostały zaobserwowane pewne wahania.
ggplot(clean_df, aes(length)) + geom_density(fill="blue")
Zmienna długości śledzia ma rozkład zbliżony do normalnego, z pewnymi perturbacjami wprowadzającymi nieregularny kształt wykresu gęstości.
Sprawdzona zostanie korelacja pomiędzy zmiennymi z wykorzystaniem współczynnika korelacji Pearsona.
corr <- round(cor(clean_df[,c(-1, -15)]),1)
ggcorrplot(corr, hc.order = TRUE,type = "lower",
lab = TRUE)
Wyraźnie najwyższa korelacja występuje pomiędzy dwoma parami zmiennych:
Z punktu widzenia zadania najważniejsza jest jednak korelacja ze zmienną length, reprezentującą długość śledzia.Najniższą ujemną korelację ma ona z temperaturą przy powierzchni wody sst (-0.4) i oscylacjami północnoatlantyckimy nao. (-0.3). Najwyższą dodatnią korelację z długością śledzia obserwujemy dla natężenia połowów w regionie (fbar) oraz dla dostępności planktonu chel1 i lcop1.
Czas w zbiorze danych nie został jawnie podany. Wiadomo jednak, że kolumna xmonth określa miesiąc, w którym dokonano pomiaru, oraz że dane w tabeli są ustawione chronologicznie latami. Z tego względu możliwe jest obliczenie średniej długości śledzia w danym miesiącu, uszeregowanie miesięcy w danym roku chronologicznie i ukazanie na wykresie długości śledzi na przestrzeni miesięcy. Oś x reprezentuje kolejne klastry miesięcy, w ktorych dokonywano pomiarów.
require(data.table)
dt_monthly <- as.data.table(clean_df)[, mean(length), by=.(xmonth, rleid(xmonth))]
df_monthly <- as.data.frame(dt_monthly)
df_monthly[,"year"] <- NA
vec <- vector()
empty_vec <- vector()
year <- 0
years <- vector()
for(i in 1:nrow(df_monthly)) {
row <- df_monthly[i,]
if(row["xmonth"] %in% vec) {
year <- year + 1
vec[vec %in% empty_vec]
}
years <- c(years, year)
vec <- c(vec, row["xmonth"])
}
df_monthly["year"] <- years
df_monthly <- df_monthly %>% group_by(year)
df_monthly <- df_monthly[order(df_monthly$year,df_monthly$xmonth),]
df_monthly["ID"] <- seq.int(nrow(df_monthly))
df_monthly
## # A tibble: 3,811 x 5
## # Groups: year [3,800]
## xmonth rleid V1 year ID
## <dbl> <int> <dbl> <dbl> <int>
## 1 3 4 22.5 0 1
## 2 5 3 25.5 0 2
## 3 6 2 23.1 0 3
## 4 7 1 23.4 0 4
## 5 6 5 23.5 1 5
## 6 3 6 22.8 2 6
## 7 4 8 22.3 3 7
## 8 6 7 23 3 8
## 9 6 9 23.2 4 9
## 10 5 10 23.7 5 10
## # ... with 3,801 more rows
ggplot(df_monthly, aes(ID, V1)) +
geom_point() +
geom_smooth() +
ggtitle("Wykres zmiany długości śledzia na przestrzeni czasu") +
ylab("Długość śledzia") +
xlab("Miesiąc pomiarowy")
Wyraźnie widzimy, że początkowo długość śledzi wzrastała, aż do osiągnięcia maksymalnej wartości, a następnie przez większość lat konsekwencje malała.
Szczegółowo długość śledzia można przeanalizować na wykresie interaktywnym wygenerowanym z użyciem Plotly:
W pierwszej kolejności zbiór danych zostanie podzielony na zbiór uczący, zawierający 80% wszystkich przykładów, oraz zbiór testowy, zawierający 20% rekordów. Dane nie zostały znormalizowane, gdyż normalizacja powodowała pogorszenie miar oceny RMSE i R2
row.number <- sample(1:nrow(clean_df[,-1]), 0.8*nrow(clean_df[,-1]))
train = clean_df[row.number,-1]
test = clean_df[-row.number,-1]
dim(train)
## [1] 42065 15
dim(test)
## [1] 10517 15
Funkcja wykorzystywana do ewaluacji modelI
test_model <- function(model, test) {
r_squared <- function (x, y) cor(x, y) ^ 2
pred <- predict(model, newdata = test)
rmse <- RMSE(pred, test$length)
r_2 <- r_squared(pred, test$length)
c(RMSE = rmse, R2=r_2)
}
Przetestowane zostaną trzy regresory:
model1 = lm(length ~., data=train)
summary(model1)
##
## Call:
## lm(formula = length ~ ., data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.1570 -0.9396 -0.0009 0.9374 6.9651
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.624e+01 7.128e+00 7.890 3.09e-15 ***
## cfin1 1.276e-01 7.583e-03 16.823 < 2e-16 ***
## cfin2 2.256e-02 3.401e-03 6.633 3.33e-11 ***
## chel1 -2.500e-03 1.420e-03 -1.761 0.078212 .
## chel2 -2.746e-03 1.943e-03 -1.414 0.157515
## lcop1 1.374e-02 1.325e-03 10.368 < 2e-16 ***
## lcop2 9.435e-03 1.594e-03 5.917 3.30e-09 ***
## fbar 6.240e+00 9.320e-02 66.946 < 2e-16 ***
## recr -3.638e-07 3.120e-08 -11.662 < 2e-16 ***
## cumf -1.038e+01 1.761e-01 -58.943 < 2e-16 ***
## totaln -5.407e-07 5.775e-08 -9.362 < 2e-16 ***
## sst -1.237e+00 2.232e-02 -55.395 < 2e-16 ***
## sal -3.803e-01 2.017e-01 -1.885 0.059414 .
## xmonth 8.276e-03 2.439e-03 3.393 0.000691 ***
## nao 2.342e-02 4.580e-03 5.112 3.20e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.372 on 42050 degrees of freedom
## Multiple R-squared: 0.3147, Adjusted R-squared: 0.3145
## F-statistic: 1380 on 14 and 42050 DF, p-value: < 2.2e-16
Model ten ma niski współczynnik determinacji R2, co znaczy, że nie jest dobrym modelem, zostanie więc podjęta próba wzmocnienia regresora poprzez podniesienie zmiennych do kwadratu.
model2 = lm(log(length)~cfin1+cfin2+chel1+chel2+lcop1+lcop2+fbar+recr+
cumf+totaln+sst+sal+nao+xmonth+ I(cfin1^2)+ I(cfin2^2)+I(chel1^2)+ I(chel2^2)+ I(lcop1^2)+
I(lcop2^2)+ I(fbar^2)+ I(recr^2)+ I(cumf^2)+ I(totaln^2 + I(sst^2)+I(sal^2)+I(nao^2)+I(xmonth^2)), data=train)
summary(model2)
##
## Call:
## lm(formula = log(length) ~ cfin1 + cfin2 + chel1 + chel2 + lcop1 +
## lcop2 + fbar + recr + cumf + totaln + sst + sal + nao + xmonth +
## I(cfin1^2) + I(cfin2^2) + I(chel1^2) + I(chel2^2) + I(lcop1^2) +
## I(lcop2^2) + I(fbar^2) + I(recr^2) + I(cumf^2) + I(totaln^2 +
## I(sst^2) + I(sal^2) + I(nao^2) + I(xmonth^2)), data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.27443 -0.03504 0.00182 0.03745 0.24635
##
## Coefficients:
## Estimate
## (Intercept) 8.422e+00
## cfin1 1.021e-02
## cfin2 -1.612e-03
## chel1 -2.180e-04
## chel2 4.015e-03
## lcop1 8.906e-04
## lcop2 -1.809e-03
## fbar 4.803e-01
## recr 2.488e-08
## cumf -5.420e-01
## totaln -7.151e-08
## sst -4.153e-02
## sal -1.305e-01
## nao -8.297e-04
## xmonth 2.457e-04
## I(cfin1^2) -2.361e-04
## I(cfin2^2) 4.962e-05
## I(chel1^2) 3.146e-06
## I(chel2^2) -7.834e-05
## I(lcop1^2) -8.126e-06
## I(lcop2^2) 3.024e-05
## I(fbar^2) -2.684e-01
## I(recr^2) -2.356e-14
## I(cumf^2) 1.360e-01
## I(totaln^2 + I(sst^2) + I(sal^2) + I(nao^2) + I(xmonth^2)) 3.475e-14
## Std. Error
## (Intercept) 3.281e-01
## cfin1 4.278e-04
## cfin2 3.206e-04
## chel1 1.332e-04
## chel2 2.107e-04
## lcop1 1.177e-04
## lcop2 1.854e-04
## fbar 1.143e-02
## recr 5.089e-09
## cumf 2.559e-02
## totaln 8.486e-09
## sst 1.017e-03
## sal 9.315e-03
## nao 2.115e-04
## xmonth 9.593e-05
## I(cfin1^2) 2.333e-05
## I(cfin2^2) 1.731e-05
## I(chel1^2) 2.172e-06
## I(chel2^2) 3.695e-06
## I(lcop1^2) 2.073e-06
## I(lcop2^2) 2.489e-06
## I(fbar^2) 1.261e-02
## I(recr^2) 3.581e-15
## I(cumf^2) 4.667e-02
## I(totaln^2 + I(sst^2) + I(sal^2) + I(nao^2) + I(xmonth^2)) 7.704e-15
## t value
## (Intercept) 25.665
## cfin1 23.858
## cfin2 -5.029
## chel1 -1.637
## chel2 19.054
## lcop1 7.564
## lcop2 -9.761
## fbar 42.010
## recr 4.889
## cumf -21.181
## totaln -8.427
## sst -40.845
## sal -14.008
## nao -3.924
## xmonth 2.562
## I(cfin1^2) -10.119
## I(cfin2^2) 2.866
## I(chel1^2) 1.448
## I(chel2^2) -21.203
## I(lcop1^2) -3.920
## I(lcop2^2) 12.151
## I(fbar^2) -21.286
## I(recr^2) -6.579
## I(cumf^2) 2.914
## I(totaln^2 + I(sst^2) + I(sal^2) + I(nao^2) + I(xmonth^2)) 4.510
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## cfin1 < 2e-16 ***
## cfin2 4.96e-07 ***
## chel1 0.10158
## chel2 < 2e-16 ***
## lcop1 3.98e-14 ***
## lcop2 < 2e-16 ***
## fbar < 2e-16 ***
## recr 1.02e-06 ***
## cumf < 2e-16 ***
## totaln < 2e-16 ***
## sst < 2e-16 ***
## sal < 2e-16 ***
## nao 8.72e-05 ***
## xmonth 0.01042 *
## I(cfin1^2) < 2e-16 ***
## I(cfin2^2) 0.00416 **
## I(chel1^2) 0.14765
## I(chel2^2) < 2e-16 ***
## I(lcop1^2) 8.86e-05 ***
## I(lcop2^2) < 2e-16 ***
## I(fbar^2) < 2e-16 ***
## I(recr^2) 4.80e-11 ***
## I(cumf^2) 0.00357 **
## I(totaln^2 + I(sst^2) + I(sal^2) + I(nao^2) + I(xmonth^2)) 6.50e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05371 on 42040 degrees of freedom
## Multiple R-squared: 0.338, Adjusted R-squared: 0.3376
## F-statistic: 894.2 on 24 and 42040 DF, p-value: < 2.2e-16
Dopasowanie zwiększyło się, ale nadal nie jest zadowalajace. Ponadto w zbiorze znajduje się wiele zmiennych o wysokich wartościach P-Value, jednak nie zostaną one usunięte, ponieważ nie przynosi to poprawy wyniku. Podjęta zostanie jednak próba zwalidowania modelu na zbiorze testowym w celu późniejszego porównania z innymi modelami:
pred1 <- predict(model2, newdata = test)
rmse <- sqrt(sum((exp(pred1) - test$length)^2)/length(test$length))
c(RMSE = rmse, R2=summary(model2)$r.squared)
## RMSE R2
## 1.3224301 0.3379539
Błąd średniokwadratorwy zwrócony przez predyktor jest znaczny w porównaniu z przedziałem wartości, jakie osiagała historycznie długość śledzia.
trControl <- trainControl(method = "cv",
number = 10,
search = "grid")
rf_default <- train(length~.,
data = train,
method = "rf",
preProc = c("center", "scale"),
ntree=10,
importance = TRUE,
metric = "RMSE",
trControl = trControl)
print(rf_default)
## Random Forest
##
## 42065 samples
## 14 predictor
##
## Pre-processing: centered (14), scaled (14)
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 37858, 37858, 37859, 37860, 37857, 37859, ...
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 2 1.169460 0.5017382 0.9234010
## 8 1.162713 0.5083092 0.9134648
## 14 1.170252 0.5026453 0.9197093
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 8.
test_model(rf_default, test)
## RMSE R2
## 1.1407368 0.5159856
Algorytm Random Forest osiąga trochę lepsze wyniki, jednak nadal są one dalekie od dobrego dopasowania.
knn_model <- train(
length ~., data = train, method = "knn",
trControl = trainControl("cv", number = 10),
preProcess = c("center","scale"),
tuneLength = 20
)
plot(knn_model)
Najlepszy wynik osiągnięto dla następującej liczby najlepszych sąsiadów:
knn_model$bestTune
## k
## 3 9
Wynik ewaluacji na zbiorze testowym:
test_model(knn_model, test)
## RMSE R2
## 1.1316467 0.5231074
Wynik jest lepszy od tych uzyskanych przy użyciu poprzednich regresorów - random forest i regresji liniowej.
Widzimy wyraźnie, że dla algorytmu K najbliższych sąsiadów największe znaczenie miała wartość temperatury przy powierzchni wody. Duże znaczenie miały też oscylacje północnoatlantyckie, natężenie połowów w regionie oraz dostępność planktonu lcop1 i chel1.
ggplot(varImp(knn_model))